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I . INTRODUCTION 


The  quality  of  tomographic  images  reconstructed  from  sparse  data  sets 
usually  suffers  from  the  presence  of  artifacts,  fuzziness  and  a  loss  of 
resolution.  Sparse  data  problems  arise  typically  in  situations  viiere  only  a  few 
views  are  available,  viiere  there  is  an  angular  restriction  in  taking  the  data, 
or  viiere  the  object  of  interest  is  partly  obstructed. 

Previous  attempts  at  sparse  data  reconstructions1  have  net  with  varying 
degrees  of  success.  Here  ve  present  a  new  algorithm  based  on  the  maximum 
2 

entropy  formalism  and  formulated  as  a  constrained  optimization  problem  which 
is  salved  by  a  finite  element  method.  The  results  show  encouraging  i irprovement 
over  conventional  sparse  data  reconstructions.  In  particular,  the  POTT 

2 

algorithm  of  Minerbo  requires  that  the  input  density  data  "falls-off  to  zero 
tovards  the  edges  of  the  target  grid.  The  present  method  lias  no  such 
restriction  on  the  input  data.  Of  course,  in  practice,  the  PENT  restriction 
means  that  the  x-ray  sources  and  detector  arrays  must  be  pulled  back,  away  from 
the  object  of  interest,  so  that  the  object  is  surrounded  by  some  medium  of 
known  density.  This  known  background  density  may  then  be  subtracted  from  the 
actual  measured  data,  resulting  in  modified  data  indicating  zero  density  near 
the  edges.  Naturally,  there  is  a  physical  restriction  resulting  from  this 
mathematical  "edge  condition”.  That  is,  given  a  fixed  scanning  installation, 
the  PENT  algorithm  may  only  be  used  to  reconstruct  density  profiles  of  objects 
whose  size  is  somewhat  smaller  than  that  which  the  machinery  itself  would 
otherwise  permit.  This  restriction  becomes  particularly  severe,  then,  in  the 
case  of  industrial  tonography,  where  the  objects  of  interest  are  already  rather 
large . 

Computational  experience  at  this  laboratory  and  elsewhere  has  demonstrated 
that  in  cases  where  data  is  sparse,  i.e.  where  fever  than  18©  views  are 
available  or  where  the  number  of  detectors  per  view  are  less  than  100,  filtered 
backprojection,  the  algorithm  of  choice  for  modern  CAT  scanners,  produces  an 
image  inferior  in  quality  to  that  obtainable  from  maximum  entropy.  The 
explanation  seems  to  revolve  around  the  built-in  smoothing  of  the  entropy 
approach  vhich  vibs  developed  on  ideas  based  on  Shannon's  information  theory. 

In  this  study  the  maximum  entropy  method  is  reformulated  and  coupled  with 
a  finite  element  approach.  This  new  algorithm  is  described  in  the  second 
chapter  followed  by  sev»eral  examples.  An  example  of  an  object  with  nonzero 
density  near  the  edge  of  the  target  grid  has  been  included.  We  demonstrate  both 
the  PENT  reconstruction  of  this  image  (which  is  thoroughly  unrecognizable )  and 
the  reconstruction  by  our  new;  algorithm,  using  the  same  measured  input  data 
The  concluding  chapter  discusses  the  results.  The  Appendix  lists  the  computer 
program  which  was  deve loped  for  this  work. 


1A.K.  Louis,  "Approximation  of  the  Radon  Transform  from  Samples  in  Limited 
Range."  in  Lecture  Notes  in  Itdical  Informatics,  Vol.  8,  pp. 127-139,  1981. 
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G.  Minerbo,  "A  Maximum  Entropy  Algorithm  for  Reconstructing  a  Source  from 
Projection  Data."  Comp.  Graph.  Irmage  Proc.  Vol.  10,  pp. 48-68,  1979. 


II.  THE  ALGORITHM 


A.  Background 

Confuted  tomography  enables  the  determination  of  density  cross-sections  of 
slices  of  an  object  in  the  plane  of  a  radiation  source  from  the  recorded 

3 

absorption  lewis  of  the  transmitted  radiation.  From  the  early  work  of  Radon 
we  know  that  the  problem  can  be  fornulated  mathematically  as  an  inverse  problem 
and  a  unique  solution,  i.e.  a  reconstructed  image  is  guaranteed  when  data  from 
an  infinite  number  of  views  is  available.  If  only  a  finite  number  of  view  are 
available,  mathematically  the  problem  becomes  ill  posed  and  small  changes,  i.e. 
discrepancies,  are  anplified  leading  to  inaccurate  image  representation. 

Techniques  haw  been  dew  loped,  however,  to  overcone  this  problem  and  to 
obtain  excellent  images.  Of  course,  image  quality  also  inproves  with  an 
increase  in  the  number  of  views.  In  this  study  we  look  at  the  case  of  21  views 
with  at  least  50  detectors  per  view.  The  problem  of  well  posedness,  existerce, 

4 

uniqueness  of  our  formulation  are  dealt  with  in  an  acconpanying  report  .  Here 
we  present  a  new  algorithm  for  the  reconstruction  of  the  image  which  for  sparse 
data  sets  yields  fewer  artifacts  than  conventional  approaches. 


B.  General  Formulation 

Central  to  the  idea  of  tomographic  reconstruction,  given  a  radiation 

source,  a  detection  system  and  an  object  placed  in  the  path  of  the  radiation, 

(see  Fig.  la)  is  the  fact  that  data  is  obtained  as  a  set  of  integrals.  Thus, 

let  f (x,y)  be  the  x-ray  attenuation  at  a  point  (x,y)  in  the  plane.  Ifeasured 

data,  G.  ,  is  available  in  the  form 
jm' 


G.  =  S.  (f)  =  M1"1*  I  f(s  cos  9.  -  t  sin  9., 

Jm  1  Jg  1  J  J 

jm 

m  =  1 , . . . m( j ) ;  j  =  1, 


s  sin  9  *  t  cos  9  .)  dtds 

J  j' 


. . J,  ( 1 ) 


where  J  is  the  number  of  projections  and  m(j)  the  number  of  detectors  for  the 
j-th  view  and  9  is  the  scan  angle. 

2 

The  object  of  interest  lies  in  a  finite  region  2)  contained  in  R  .The 
entropy  of  the  image  can  be  defined  as 


J.  Radon,  "Ueber  die  Best  irnrung  von  Funktionen  durch  lhre  Integralwerte  laengs 
gewisser  Ifennigfalt igkeiten. “  Ber.  Verh.  Saechs.  Akad.  Wiss.  Leipzig,  Math. 
Phys.  Kl.  Vol.  69,  pp.  262-277,  1917. 

*R.T.  Smith,  C.K.  Zoltani,  "An  Application  of  the  Finite  Element  ffethod  to 
Maximum  Entropy  Tomographic  Image  Reconstruction."  ERL  Report  (in 
preparation) . 


where  the  constraints  set  L  is  defined  as 


f  =  {  f  5  2  ):  a  <  f  <  b,  and  f  =  0  in  fR2  \  2)}  (5) 


vjhere  we  usually  require  a  >  0  and  b  <  co, 

2 

From  the  theory  of  penalty  functions,  we  knovi  that  if  we  take  L  =  i.  (  3)  ) , 
then  as  y->  to,  the  solution  of  the  unconstrained  minimization  problem  converges 

2 

to  the  solution  of  the  equality  constrained  problem  solved  by  Miner bo  . 

By  taking  a  sufficiently  large  value  of  the  penalty  paraiieter  Y,  the  residual 

error,  >  {G  -  S  (f))^  can  be  made  sufficiently  small  to  achieve  the  desired 

’  Z  J1"  J™ 

fidelity  in  the  reconstructed  image.  In  practice,  of  course,  the  values  of  G 

jm 

are  degraded  by  noise  and  we  wish  only  to  obtain  a  total  residual  error  within 
some  tolerance  determined  by  the  known  accuracy  of  the  measurements.  Also,  we 
have  included  a  priori  information  in  the  problem  formulation,  by  choosing  a 
and  b  so  that  the  attenuation  lies  within  sone  known  physical  limits. 
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Figure  1.  Scanning  Geonetry 
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The  solution  of  (4)  can  be  approx imated  by  solving  the  problem  in  a  finite 

2  Yi  jj 

dimensional  subspace  of  &  (  3D  ),  S  \  Let  {4^(K)  «  *  •  •  »4'n(x)}  be  a  basis  for  S 


for  any  c  =  (c^  ,c2, . .  ,0^)  in  Pt^let 


fc(*)  =  ^  Ci‘MM>-  (6) 

i=l 

It  is  shown  in  [4]  that,  for  a  certain  class  of  approximating  subspaces 

(including  the  one  used  in  the  present  work),  the  solution  of  (4)  out  of  Sh 
converges  to  the  solution  of  the  infinite  dimensional  optimisation  problem,  as 
h  tends  to  zero.  Precise  estimates  of  the  deviation  of  the  reconstructed  image 
from  that  of  the  true  phantom  are  out  of  the  question  at  this  tint;,  as  an 
analytic  relationship  between  the  maximum  out ropy  solution  and  the  true  image 
is  unknown. 

The  finite  dimensional  constrained  optimization  problem  then,  is  to 

determine  f  which  minimizes 
c 

V  V  V 

E(£  )  =  -  i?(  >  c  4*.  )  +  y  >  (G .  -  S  (  >0.4*  ))*  (?) 

'  c'  L  iTl'  L  v  J**1  tnr£  iti'/  '  1 


subject  to 


<  ^  ci4’i(x)  £  b,  for  ill  x  £  B  and  a  >  0. 


To  solve  this  nonlinear,  convex  programming  problem,  a  uniform  mesh  is 
superimposed  on  the  region  2)  with  step  size  h  in  both  coordinate  directions,  f+ 
product  of  piecewise  linear  functions,  4*  (x»y), 


4‘ij(x,y)  =  4»i(x)4,j(y ) »  1  =  1*  -  -  -m1 

j  =  1 ,  .  .  . 


is  used  as  a  basis  for  S  .  Here,  denotes  the  usual  linear  finite  e 


lenent 


basis  function 


0,  if  x  <  x  ,  or  x  :  x . , , 
l-l  l+l 

4*  (x)  -  (x  -  x  .  )/h  if  x  .  x  x 

fi'  ’  '  i-l'  i-l  j 

(x.^-xj/hifx  i  x  i  v.  , 

'  l+l  *  i  i+l 


vjhere  k  . 

1 


ih. 


For  convenience,  we  write 

=  ^(xH’jfy).  k  =  1  +  (j  -  l)m,  (10) 

then 

M 

fc(x>y)  =  ^  ciA(*!y) 

k=i 

.here  II  -  rr|m2  ’ 

idling  the  compact  support  of  the  s,  it  can  be  shown  that 

a  <  fc(x,y)  <  b  for  all  (x,y)  f  39 

if  and  only  if  a  <  c^<  b,  for  all  k  =  1,...M. 

The  optimization  problem  reduces  to  one  with  linear  inequality 
constraints,  i.e.  determine  c  which  minimizes 

E(f  )  =  -r?(F  )  +  V  -  S  .  (f  ))2  (11) 

'  c'  '  c1  jm  jm'  c' '  1 

J,m 

subject  to  a  <  Cj,  <  b,  k  =  1,...M. 

T 

The  problem  then  is  to  calculate  the  vector  c  =  (c^,...^)  for  gi<»en  valiES  of 
V,  a  and  b. 

For  the  optimization  step  of  the  algorithm  we  used  MINOS5,  a  FORTRAN  code 
developed  for  constrained  optimization  problems.  The  user  supplies  the  function 
to  be  optimized,  the  constraints  and  some  ancillary  files. 

The  calculation  proceeds  as  follow:  a  preprocessor  program  (see  Fig.  2a) 
reads  the  number  of  nodes  in  the  mesh  geometry,  the  number  of  projection  angles 

B.A.  Murtagh,  M.A.  Saunders,  "MINOS  5.0  User’s  Guide."  Systems  Optimization 
Laboratory,  Department  of  Operations  Research,  Technical  Report  SCI.  03  20, 
Stanford  University,  Stanford,  CA,  December,  1983. 
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and  their  values  (the  ANG.DAT  file),  and  the  nunt>er  of  detectors  for  each 
projection  angle  (view).  The  preprocessor  creates  a  weights  file  SJM.DAT,  which 
indicates  the  contribution  of  each  particular  finite  element  (geometry  cell)  to 
each  detector  at  each  angle. 

The  main  program  is  the  optimizing  routine  MINOS  (see  Fig.  2b).  ft  first 
guess  at  the  coefficients  c^  is  read  into  MIN06.  On  the  first  pass  through,  the 

subroutine  FUNOBJ  reads  the  measurements  file  GJM.DAT,  which  contains  the 
attenuation  values  indicating  hov;  much  of  the  incident  beam  was  absorbed  by  its 
traverse  through  the  target.  These  measurements  values  are  scaled  such  that  the 
maximum  attenuation  is  set  to  1.0.  On  each  successive  iteration,  FUNOBJ 
accesses  the  weights  file  SJff.DAT  in  its  calculation  of  the  entropy  EXT  for  the 
present  values  of  the  coefficients  c^  and  the  terms 


M 

EN(k)  =  J  J  <^(K,y)  ln[^  ck<^(x,y)]  dxdy  (12) 

SO  k=i 

which  are  used  in  computing  the  gradient  of  the  object  function.  For  the 
quadrature  of  the  equation  (12)  the  adapt ive  grid  routine  QUANC8&  was  used. 

With  the  value  of  EXT  and  EN(k),  FUNOBJ  will  compute  the  values  of  the  object 
function  F(c)  and  the  gradient  function  G  where 


M  M 

F(c)  =  -Ol^c^lK.y))  *  [Gjm- 

k=l  j,m  k= 1 

and 

M  M 

g(c)  =  v  f(c)  =  ^>ek{h2  +  J  J  lnt^  ck^(K»y)3  dKrfy 

k=l  E*  k=l 

M 

j,m  k=l 


^G.E.  Forsythe,  M.  A.  Malcolm,  C.B.  Ifoler,  Conputer  Methods  for  Mathemat ical 
Computations,  Prentice  Hall,  1977. 
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where  efc  is  the  k-th  unit  vector,  i.e.  efc  =  (0,®, . . .  1, . .  .®)T  are!  E*1  is  the 

support  of  the  k-th  basis  function,  4^(x,y).  The  values  of  F  and  G  are  then 

passed  back  to  HINOS  for  the  next  step  in  the  optimization.  The  value  of  the 
penalty  parameter  v  is  problem  dependent  and  some  experimentation  is  reeded  to 

7 

determine  the  best  value. 

The  solution  can  be  made  arbitrarily  close  to  the  optimal  solution  by  choosing 
a  large  enough  value  for  y.  An  analytical  nethod  for  determining  tte  optimal 
penalty  is  not  available.  Also,  some  caution  needs  to  be  exercised  in  picking 
the  value  to  avoid  ill-conditioning  or  slow  convergence.  In  this  problem  the 


penalty  parameter  was  increased  gradually  and  it  was  found  that  y  =  75  yielded 
the  best  result,  i.e.  the  smallest  error  in  the  reconstruction. 


III.  IPFLETENTATION  OF  TOE  TECHNIQUE 


A.  ffathenat ical  Phantoms 

The  computer  code  GOLEM8  was  used  to  generate  the  x-ray  absorption  data 
from  mathematical  phantoms  for  input  to  the  reconstruct  ion  algorithms.  For  this 
study,  three  examples  were  generated,  assuming  a  parallel  mode  of  scanning, 
using  a  technique  described  in  Ref.  8.  Physically,  the  target  objects  are 
placed  between  a  source  of  monochromatic  parallel  beam  x-rays  and  a  straight 
line  of  detectors  perpendicular  to  the  x-ray  transmission.  For  each  of  the 
examples,  the  phantoms  were  set  within  a  square  target  grid  of  30  unit  cells  on 
each  side.  The  x-ray  source  was  placed  69  units  from  the  center  of  the  target 
grid.  A  line  of  25  detectors  was  placed  16  units  from  the  center  of  the  target 
grid,  opposite  the  x-ray  source.  Data  was  produced  from  five  projection  angles 
equally  spaced  around  the  target.  A  square  reconstruction  grid  of  30  unit  cells 
on  each  side  spans  the  width  of  the  25  detector  bins. 

The  first  phantom  consisted  of  a  tube  6  units  in  diame+^r  ,  with  unit 
density  on  a  0.1  density  background,  whose  center  was  placed  on  the  axis  of  the 
field  of  scan.  For  the  second  sample  problem,  three  identical,  solid  cylinders 
were  placed  at  arbitrary  positions  within  the  target  area.  In  the  third 
problem,  a  solid  cylinder  was  placed  inside  a  thick-walled  hollow  cylinder, 
with  both  cylinders  centered  at  the  target  origin.  For  all  the  examples,  the 
cylinders  were  given  a  density  of  1  while  the  remainder  of  the  target  grid  was 
assigned  a  background  density  of  0.1.  A  cross-sectional  view  of  the  sample 
problems  is  given  in  Figure  1. 
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B.  Results 


The  sanple  phantoms  were  reconstructed  using  our  new  algorithm  and  then 
conpared  with  the  results  of  the  currently  preferred  sparse  data  technique 

PENT2.  Figure  4a  shows  the  reconstructed  density  plot  of  the  first  exanple 
using  the  new  technique,  while  the  adjoining  Figure  4b  illustrates  the  PENT 
results  for  the  sane  problem.  These  figures  are  both  three-dimensional  grid 
representations  of  the  reconstructed  density  (with  hidden  lines  removed).  The 
31x31  grid  used  for  the  reconstruction  (the  actual  finite  element  grid 
for  Figure  4a)  is  shown  in  the  middle  of  the  51x51  grid,  with  points  outside 
the  31x31  grid  being  given  density  zero.  One  should  note  the  high  degree  of 
correspondence  between  the  reconstructions  of  Figures  4a  and  4b. 

Figure  5a  gives  the  reconstructed  density  plot  of  the  second  exanple.  The 
three  solid  cylinders  of  unit  density  stand  out  vividly  against  an  undulating 
background  of  density  0.1.  The  five  projection  angles  can  be  picked  out  as  the 
'ridges''  or  "hunps"  which  stretch  through  the  object  space  in  Figure  5a.  One 
should  notice  from  Figure  3b  that  one  of  the  three  cylinders  is  rather  close  to 
the  edge  of  the  target  grid.  The  PENT  algorithm  has  difficulty  w/ith  nonzero 
density  near  edges.  In  Figure  5b  one  can  see  the  PENT  reconstruction  of  the 
second  sample  problem  using  the  same  projection  data  as  that  used  in  Figure  5a. 
Of  course,  the  reconstruction  is  completely  unrecognizable.  In  Figure  5c  we 
demonstrate  the  results  of  a  MENT  run  in  which  the  data  from  the  three  cylinder 
geometry  case  is  reconstructed  in  a  larger  area  (a  51  x  51  grid).  Pbving  the 
cylinders  away  from  the  edge  of  the  reconstruction  grid  allows  the  PEKT  code  to 
smooth  its  results  to  zero  at  the  boundaries  of  the  reconstruction  grid.  In 
particular,  the  PEKT  algorithm  artificially  forces  the  density  to  zero  near  the 
edges  of  the  target  grid,  while  the  new  method  correctly  shows  a  nonzero 
Jensity. 

The  results  for  the  third  sanple  problem  are  given  in  Figures  6  and  7. 
Figure  6  gives  density  surface  plots,  comparing  the  nevi  algorithm  vrith  PEKT. 
Finally  in  Figure  7,  a  vertical  slice  through  tbs  center  of  the  density  plot  is 
shown  for  both  methods.  In  this  case,  note  that  there  is  very  good  agreement 
between  the  two  net  hods.  However,  as  with  the  second  exanple,  PEKT  lias  forced 
the  density  to  zero  near  the  edges. 

The  new  algorithm  was  run  on  the  CRAY  XP1P/12  computer  at  the  Naval 
Research  laboratory  in  Uashington,  DC.  The  results  shown  were  obtained  using 
approximately  100  iterations  of  the  P1IN06  optimization  code  and  required 
approximately  10  minutes  of  cpu  time.  The  PENT  algorithm  was  run  on  the  CDC 
CYBER  70/76  conputer  at  the  Ballistic  Research  Laboratory,  Aberdeen  Proving 
Ground,  ft).  These  results  required  only  a  few  seconds  of  cpu  time. 

Although  the  new  method  requires  substantially  more  conputer  time  and 
memory  than  the  PENT  algoritlim,  we  believe  that  it  still  represents  a 
significant  advance.  Very  sparse  data  is  ail  tliat  is  available  in  many 
industrial  applications.  One  can  enploy  no  mathematical  trickery  to  change  this 
fact.  In  some  of  these  cases,  where  the  data  fails  to  fall  off  to  zero  near  the 
edges,  the  currently  preferred  method  PENT  will  not  produce  an  even  remotely 
recognizable  reconstruction.  Our  new  algorithm  has  removed  this  "edge 
condition"  restriction  and  altliough  a  costly  process  (in  terms  of  conputer 
time),  it  wrill  yield  a  reasonable  reconstruction.  We  believe  that  the  conputer 


tine  and  memory  requirements  nay  be  reduced  by  developing  more  efficient  codes 
for  the  computation  of  the  cost  functional  and  gradient  terms  in  (13)  and  (14) 
and  possibly  by  the  developnent  of  a  customized  optimization  routine. 


b.  Sanple  Problem  Two 

c.  Sanple 

Problem  Three 
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Figure  3. a.  Object  Geometry  Sanple  Problem  One 

b.  Object  Geometry  Sample  Problem  Two 

c.  Object  Geometry  Sample  Problem  Three 


Figure  7.  a.  Sanple  Problem  Number  Two:  Center  Slice 
b.  Sanple  Problem  Nunfeer  Two:  Center  Slice. 
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SUBROUTINE  FUNOBJ(MODE,NN,C,F,G,NSTATE,NPROB, Z.NWCORE) 

THIS  SUBPROGRAM  WILL  COMPUTE  THE  OBJECT  FUNCTION  (F) 
AND  ITS  GRADIENT  (G(K) )  CORRESPONDING  TO  THE  CURRENT 
VALUES  OF  THE  COEFFICiefTS,  C(K)  COMPUTED  BY  MINOS. 

THESE  VALUES  ARE  THEN  PASSED  BACK  TO  MINOS.  TO 
CONTINUE  THE  OPTIMIZATION  PROCESS. 


DIMENSION  C(NN)  , EN(289)  ,C(HN)  ,GJ!1(231  j ,  bJfflAT'SM&a) 

LOGICAL  FIRST 
DATA  FIRST/. TRUE. / 

URITE(6,37) 

37  FORMAT}  ’  THE  COEFFICIENTS  C(K)') 

URI1E(6,35)  (C( INDEX) ,INDEK=1, 289) 

35  FORMAT  ( 17F7. 3) 

IF  (FIRST)  THEN 
UBrnE}6,22) 

22  FORMAT}  ’  PERFORMING  SET -IJP  OPERATIONS  IN  FUNOBJ’) 

OPEN  ( 4 ,  F ILE-  *  SJMDAT-‘ ,  ACCESS=  *  SEQUENTIAL  * ,  FORM-  ’  UNFOPJii  I7ED ' 

READ  (4)  SJM3AT 

OPEN  (7,FIL£=’GJMDAT‘  ) 

REWI ND ( 7 ) 

N-6 
It  11 
JA=5 

IND=(2*N+l)**2 

H=«./N 

N  =  NUMBER  OF  GRID  POINTS  ON  EITHER  SIDE  OF  THE  ORIGIN 
M  =  NUMBER  OF  PROJECTION  STRIPS 

JA  =  THE  NUMBER  OF  PROJECTION  ANGLES,  I.F.  #  OF  UIF.VF 
IND  =  NUMBER  OF  FINITE  ELEMENTS 
H  =  MESH  WIDTH 

THE  GJM(N)  ’S  ARE  THE  MEASUREMENTS  OF  THE  K-RnY  INTENSITY. 
THESE  ARE  STORED  IN  THE  FILE  GJM.DAT. 

DO  492  JJA=i , JA 
DO  492  MA=1 ,M 

492  READ(7,493)  GJM( (JJA-1 )*M+MA) 

493  FORMAT  ( El  5 . 3 ) 

CLOSE  ( 7 ,  STA'I  US  =  ’  DELETE '  ) 

WRITE} fa, 38) 

38  FORMAT} »  THE  MEASUREMENTS  FILE’ ) 

DO  33  MA=1,M 

33  WRITE(6, 495)  (GJM(  (JJA-1  )#M+MA) ,  JJA=1 ,  JA) 

495  FORMAT} 1P5E15. 3) 
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RELEER  AND  ABSERR  ARE  THE  RELATIVE  AND  ABSOLUTE 
ERRORS  ALLOWED  IN  THE  NUMERICAL  INTEGRATIONS 


RELERR=5E-6 

ABSERR=0. 

QAItIA  IS  THE  PENALTY  PARAMETER. 

FOR  CONVERGENCE  TO  THE  MAXIMUM  EKTRQPY  SOLUTION, 

UE  NERD  TO  TAKE  A  SEQUENCE  OF  VALUES  OF  GAlflA 

TENDING  TO  INFINITY,  OR  TAKE  A  SUFFICIENTLY  LARGE  VALUE 

OF  GAM1A  SO  THAT  THE  RESIDUAL,  RES  IS  SUFFICIENTLY  SMALL. 


GAHMA=50.0 
l*ITE(6,39)  GAM1A 
39  FORMAT(’  GAMMA=',F5.1) 

FIRST=. FALSE. 

URITE(6,24) 

24  FORMAT ( ’  FINISHED  UITH  SET-UP  PROCEDURES  IN  FUNOBJ') 

ENDIF 

CALL  ENIR(EN, EXT, N.RELERR, ABSERR, C) 

THE  SUBROUTINE  ENTR  COMPUTES  THE  VALUES  OF  THE  ENTROPY 
TERMS  EN(K)  AND  THE  ENTROPY  EXT  FOR  THE  CURRENT  VALUES  OF  C(K) 

RES  =49. 

DO  420  JTA=1,JA 
DO  420  Pn=l,M 
SUM=0. 

DO  410  K=1 , IND 

JJ  =  (K-1)»M*(JTA-1)*IND*M+MI 
32  FORMAT(  ’  JJ=*,I10) 

MAI W. FOR  COMPUTES  THE  VALUES  OF  SJM(K),  BUT  STORES  THESE 
VALUES  IN  SJM.DAT  IN  A  DIFFERENT  ORDER  FROM  THAT  IN  WHICH 
THEY  MUST  BE  ACCESSED  BY  FUN.  JJ  GIVES  THE  DESIRED  RECORD 
NUMBER  IN  THE  DIRECT  ACCESS  FILE  SJM.DAT. 

10  SUrfc5UM*C(K)»SJl«AT(JJ) 

20  RE5=RES+(GJM(  ( JTA-1 )*M+rtI )  -SUM)** 2 

RES  IS  THE  RESIDUAL  ERROR  IN  METING  THE  MEASUREMENTS. 

F  *  -EXT  ♦  QA!TO»RES 
WRITE(6,34)  F, EXT, RES 

34  FORMAT^  F=*,E14.5,»  EXT=' ,E14. 5, '  RES=',E14.5) 
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THIS  SECTION  COfPUTEB  THE  GRADIENT  TERMS 


DO  440  K=1 , IKD 
SU2=0. 

DO  460  JTA=1 ,  JA 
DO  460  MI=1,M 

sun=o. 

DO  450  KI=1, IND 

JJ=(KI-1 }#M+(JTA-1  )#IND#M+MI 
sun=suw+c  ( K I )  «SJM)AT  ( J  J ) 

J2=(K-1  )«M+(JTA-1  )#IND#M+MI 
SU2  =SU2+  ( GJJ1(  ( JTA-1)  *HWH )  -SUM) *SJ1CAT  ( J  2 ) 
G(K)  =-H»H*(l.+AL0G(4.))  +  EN(K)-2.*GATt1A#SU2 
RETURN 


800  PRINT  *,  *  ERROR/FUNOBJ  -  EOF  ON  UNIT  4’ 

STOP 

810  PRINT  «,  *  ERROR /FUNOBJ  -  ERR  ON  UNIT  4* 

PRINT  *,  *  IOSTAT  =  * ,  IONUM 

STOP 

END 


SUBROUTINE  ENIRfEN.ENT.N.RELERR, ABSERR.C) 

NUTERICAL  COfFUTATION  OF  ENTROPY  TERMS 

DIMENSION  C(*),EN(«) 

IND=(2*N+1)**2 

IND  GIVES  THE  NUflBER  OF  FINITE  ELEMENTS  IN  THE  GRID 

DO  50  R=1 , IND 

L=( (K-l )/(2#N+l ) )-N 
I=K-(L+N)«(2*N+1)-N-1 

I*H  AND  L#H  GIVE  THE  X  AND  Y  COORDINATES  OF  THE 
CENTROID  OF  THE  FINITE  ELEMENT,  RESPECTIVELY 


COWUTE  ENTROPY  QUADRANT-BY -QUADRANT 

INITIALIZE  ENTROPY 

EN(K)  =  0. 

ENT  =  0. 

ENT  IS  THE  ENTROPY.  EN(K)  IS  THE  CONTRIBUTION 
TO  THE  ENTROPY  FROM  THE  KTH  FINITE  ELEMENT 
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C  FIRST  QUADRANT 

C 

IF  (I.EQ.N.OR.L.EQ.N)  GOTO  60 

CALL  QUAD(N, I  ,L,C,K,  AN,RELERR,  ABSERR, 1) 

EN(K)  =  EN(K)+AN 
C 

C  SECOND  QUADRANT 

C 

60  IF  (I.EQ.-N.OR.L.EQ.N)  GOnTO  70 

CALL  QUAD(N, I,L,C,RIAN,RELERR,ABSERR,2) 

EN(K)  =  EN(K)+AN 
C 

C  THIRD  QUADRANT 

C 

70  IF  (I.E3Q.-N.OR.L.EQ.-N)  GOTO  80 

CALL  QUAD(N,  I,L,C,K,AN,RELERR,ABSERR,3) 

EN(K)  =  EN(K)+AN 
C 

C  FOURTH  QUADRANT 

C 

80  IF  (I.EQ.N.OR.L.EQ.-N)  GOTO  90 

CALL  QUAD(N,I,LIC,R,AN,RELERR,ABSERR,4) 

EN(K)  =  EN(K)+AN 
90  ENT  =  ENT+C(K)#EN(K) 

50  CONTINUE 

RETURN 
END 
C 

C  FUNCTION  SUBPROGRAM  FUR  INTEGRAND 

C 

FUNCTION  FN(X,  J,  I,K,N,C,H) 

DI FENS ION  C(w) 

IF  (J.EQ.l.OR.J.EQ.4)  GOTO  210 
F=X-(I-1.)#H 

GOTO  220 

210  F=( 1+1 . )*H-X 

C 

C  THE  VALUE  OF  J  WILL  DETERMINE  WHICH  QUADRANT  OF  THE 

C  FINITE  ELEMENT  WE  ARE  LOOKING  AT. 

C 

220  IF  (J.EQ.l)  GOTO  230 

IF  (J.EQ.2)  GOTO  240 
IF  (J.EQ.3)  GOTO  250 
GOTO  260 

AJ  =  (C(K*NM)-C(K+1))*(X-I*H)  +  (C(K+N)-C(K))*((I+1. )*H-X) 
CJ  =  C(K+i)*(X-I#H)+C(K)#((I+l. )#H-X) 


230 


240 


GOTO  309 

AJ  *  (C(R+N-1)-C(R-1))»(I*H-X)*(C(K+N)-C(R))*(X-(I-1. )*H) 
CJ  =  C(H-l)»(I«H-X)*C(R)*(X-(I-l. )»H) 

GOTO  300 

250  AJ  =  (C(R-N-1)-C(R-1))*(I«H-X)*(C(R-N)-C(R))»(X-(I-1.)*H) 

CJ  =  C(R-1)»(I«H-X)*C(R)«(X-(I-1.  )*H) 

GOTO  300 

260  AJ  =  (C(R-N+1)-C(R*1))«(X-I*H)*(C(R-N)-C(R))#((I+1.)*H-X) 

CJ  =  C(R+lj»(X-I«H)+C(E)«((I*l.)*H-X) 

300  IF  (AJ.EQ.0. )  GOTO  310 

IF  (ABS(AJ) .LT. IE-6)  GOTO  310 
BJ=AJ*CJ 

G  =  BJ*BJ«jOG(BJ)/(2.«W*AJ)-CJ»(AJ*BJ)#LjOG(CJ)/(2.«AJ*AJ) 

D  -CJ/{2.*AJ) 

GOTO  320 

310  G  =  .75  +  .  5*LOG(CJ ) 

320  FN  =  F  *  G 

RETURN 
END 
C 
C 

SUBROUTINE  QUAD(N,  I.L.C.K.AN.RELERR, ABSERR.J) 

C 

C  QUAD  SETS  UP  THE  PARAMETERS  FOR  THE  NUfERICAL  QUADRATURE 

C  QUANC8  PERFORTE  AN  8-POINT  ADAPTIVE  GRID  QUADRATURE. 

C 

DIPENSION  C(«) 

H=l./N 

IF  (J.EQ.l.OR.J.BQ.4)  GOTO  110 
A  =  (1-1. )*H 
B  =  I*H 
GOTO  115 

110  A  =  I*H 

B  =  (1*1. )«H 

115  CALL  QUANC8( A, B, ABSERR,RELERR, RESULT, ERREST , NOFUN . FLAG, 

D  J, I ,K,N,C,H) 

IF  (FLAG.NE.0. )  URITE(«,2)  FLAG 

2  FORMAT  (48H  WARNING. .RESULT  MAY  BE  UNRELIABLE.  FLAG  =  ,F6.2) 

AN  =  RESULT* (H*H/4.  )*(L0G(4.  )-LOG(H) )-3.*H*H/B. 

RETURN 

EM) 
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SUBPROGRAM  FOR  ADAPTIVE  GRID  QUADRATURE 


SUBROUTINE  QUANC8(A,B, ABSEMi,RELERRf RESULT, EEREST, NOFUN, FLAG, 
D  J , I ,K,N,C,H) 


DIMENSION  QRIGHT(31 ) ,F( 16) ,X(16) ,F5AVE(8,30) ,XSAVE(8,30) ,C(*) 

LEVMIN=i 

LEVHAX=30 

LEVOUT=6 

NOMAX=5000 

NOFIN=NOHAX-8*(LEVI1AX-LEVOUT+2<h*(LEVOUT+1  ) ) 

U0=3956./14175. 

Wl=23552./14175. 

U2=-3712./14175. 

113=41984.  / 14175. 

W4=-18160./14175. 


FLAG=0. 

RESULT-0. 

CC«11=0. 

ERREST=0. 

AREA=0. 

NOFUN=0 

IF  (A.EQ.B)  RETURN 


LEV=0 
NIIt=l 
X0=A 
X(  16)=B 
QPREV=0. 

FteFNfXC.J.I.K.N.C.H) 

ST0NE=(B-A)/16. 

X(8)=(X0+X( 16) )/2. 
X(4)=(X0+X(8))/2. 

X( 12)=(X(8)+X( 16) )/2. 
X(2)=(X0+X(4) )/2. 
X(6)=(X(4)+X(8) )/2. 
X(10)=(X(8)+X(12) )/2. 
X(14)=(X(12)+X(16))/2. 

DO  25  JT=2,16,2 

25  F( JT)=FN(X( JT) ,J,I,K,N,C,H) 

N0FUN=9 
C 


fcj  u  '■«  LI  III  tv».  4  A 


X(l)=(X0+X(2))/2. 

F(  1 )=FN(X( 1 ) , J, I ,K,N,C,H) 

DO  35  JP=3, 15,2 

X(JP)  =  (X(JP-i)+X(JP+l))/2. 

F(JP)=FN(X(JP) , J, I ,K,N,C,H) 

N0FUN=N0FUN+8 

STEP=(X(16)-X0)/16. 

QLEFT=  ( W* ( F0+F( 8 ) )  +W1* (F(  1 )  +F{  7 ) )  +W2* (F(  2 )  +F( 6 ) ) 

D  +W3*(F(3)+F(5))+W4«F(4))*STEF 

QRK3fr(LEW+l)  =  (W0*(F(B)+F(16))+Wlw(F(9)+F(15))+W2*(F(10)+F(14)) 
D  +W3#(F{  11)+F(  13)  )+W4#F(  12)  )»STEF 

QNOW=OLEFT+QRIGHT  ( LEV+ 1 ) 

(^)IFF=(^OW-QPREV 

AREA=AREft+QDIFF 


ESTERR=ABS(QDIFF)/1023. 

TOUSER=Aimi(ABSEXR,RElERR#ABS(AREA) )*  (STEP /STONE) 
IF  (LEW.LT.LEWMIN)  GOTO  50 
IF  (LEV.  GE.LEVTlflX)  GOTO  62 
IF  (NOFIJN.GT.NOFIN)  GOTO  60 
IF  (ESTERR.LE.TOLERR)  GOTO  70 


Nin=2KNin 

LEV=LEV+1 

DO  52  IT=1 ,8 

FSAVE( IT,LEV)=F( IT+8) 
XSflVE( IT, LEW) =X( IT+8 ) 

QPREV=QLEFT 
DO  55  IQ=1,8 
JQ=-IQ 

F(2*JQ+18)=F( JQ+9) 
X(2*JQ+18)=X(JQ+9) 

GOTO  30 

N0FIN=2*N0FIN 

LEVMflX=LEVOUT 

FLAG=FLAG+(B-X0)/(B-A) 

GOTO  70 

FLAG=FLAG+1 . 

RESULT = RESULT  -KJNCH4 
EXREST =ERRE5T+ESTERR 
COR1  l=CORl  1  +<?DIFF/ 1023 . 


E 


c 

72  IF  (NIM.EQ.2tt(NIM/2) )  GOTO  75 

HIIfcHIM/2 
L£V=LEV-1 
GOTO  72 

75  MIM=Mimi 

IF  (LEV.LE.0)  GOTO  80 
C 

QPREV=QRIGHT(LEV) 

X0=X( 16) 

F0=F(16) 

DO  78  IL=1 , 8 

F(2#IL)=FSAVE(IL,LEV) 

78  X(2«IL)=XSAVE(IL,LEV) 

GOTO  30 
C 
C 

80  RESULT  ^RESULT  +COR 1 1 

C 

C 

IF  (ENREST.EQ.0.  )  RETURN 
82  TETF=ABS  (RESULT) +ERREST 

IF  (TEJF.NE.ABS(  RESULT) )  RETURN 

ENREST=2.»EKRE5T 

GOTO  82 

END 
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